def a(n):
    return sum([1.0/e for e in range(1,n)])

def b(n):
    return sum([1.0/e**2 for e in range(1,n)])

def watterson(s,n):
    return float(s)/a(n)

def var_watterson(s,n,theta):
    return 1.0/(a(n)**2)*(s+theta**2*b(n)) 

theta=watterson(26,5)
print "theta=",theta
vartheta=var_watterson(26,5,theta)
print "var theta=",vartheta
